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I.  Barrodale  and  R.E.  Erickson 


Experience  with  the  maximum  entropy  method  of  spectral  analysis 


suggests  that  (a)  it  can  produce  inaccurate  frequency  estimates  of 


short  sample  sinusoidal  data,  and  (b)  it  sometimes  produces  calculated 


values  for  the  filter  coefficients  that  are  unduly  contaminated  by 


rounding  errors.  Consequently,  it?  this  report  sty  developsjan  algorithm 
for  solving  the  underlying  least-squares  problem  directly,  without 


forcing  a Toeplitz  structure  on  the  model.  This  approach  leads  to  more 


accurate  frequency  determination  for  short  sample  harmonic  processes 


and  our  algorithm  is  computationally  efficient  and  numerically  stable 


The  algorithm  can  also  be  applied  to  two  other  versions  of  the  linear 


prediction  problem.  A FORTRAN  program  is  supplied 
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1*  Introduction 

The  main  purpose  of  this  report  is  to  describe  an  algorithm,  and 
present  a FORTRAN  program,  for  solving  three  variants  of  the  so-called 
linear  prediction  problem.  Our  algorithm  calculates  an  "optimal" 
least-squares  (L S)  solution  by  solving  a sequence  of  normal  equations 
in  a numerically  stable  manner,  and  it  is  efficient  in  both  computer 
storage  and  time  requirements.  In  particular,  it  can  be  used  to  deter- 
mine the  parameters  of  the  autoregressive  (AR)  model  associated  with 
the  maximum  entropy  method  (MEM)  of  spectral  analysis,  and  this  applica- 
tion provided  our  primary  motivation  for  developing  the  program  given 
in  Appendix  E.  However,  the  program  should  also  prove  to  be  useful  in 
other  applications,  such  as  the  design  of  adaptive  signal  whiteners, 
estimating  the  parameters  of  adaptive  control  systems,  and  time  series 
modelling. 

We  also  discuss  some  techniques  for  preserving  accuracy  in  LS 
floating-point  calculations  which,  although  they  are  of  fundamental 
importance  in  calculating  satisfactory  solutions,  appear  to  have  been 
neglected  by  many  users.  In  addition,  we  briefly  review  some  alternative 
methods  to  handle  problems  for  which  our  algorithm  is  less  efficient. 

Currently,  the  most  popular  algorithm  for  determining  the  MEM 
parameters  is  due  to  Burg  (1967,  1968,  1975),  and  a perusal  of  the 
recent  literature  in  geophysics  alone  confirms  that  MEM  spectral  esti- 
mates obtained  with  this  algorithm  have  gained  wide  acceptance.  The 
computational  efficiency  of  this  algorithm  stems  from  its  use  of  the 
Levinson  (1947)  recursive  properties  of  solutions  to  sequences  of  certain 


linear  algebraic  equations  with  coefficient  matrices  of  symmetric  Toeplitz 
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form.  (A  Toeplitz  matrix  has  all  its  elements  equal  along  each  diagonal, 

i.e.  c.  , » c.  . .) 
i.j  1-3 

However,  when  processing  real  data  (e.g.  geomagnetic  micropulsa- 
tion recordings)  we  have  experienced  two  shortcomings  inherent  in  Burg's 
algorithm.  The  first  of  these  is  that  it  can  produce  inaccurate  frequency 
estimates  of  short  sample  sinusoidal  data.  This  limitation  appears  to 
have  been  first  reported  by  Chen  and  Stegen  (1974).  Erickson  (1976), 

Ulrych  and  Clayton  (1976),  and  Barber  and  Taylor  (1977)  provide  empirical 
evidence  which  suggests  that  LS  solutions  to  the  AR  model  for  MEM  lead  to 
more  accurate  frequency  determination  in  these  cases.  Burg's  algorithm 
forces  a Toeplitz  structure  on  the  matrix  of  the  system  of  equations 
which  yields  the  AR  parameters,  and  the  resulting  spectral  estimates  for 
short  sample  harmonic  processes  are  usually  inferior  to  the  spectral 
estimates  resulting  from  our  LS  algorithm. 

The  second  shortcoming  of  Burg's  algorithm  is  that  it  sometimes  pro- 
duces calculated  values  for  the  parameters  that  are  unduly  contaminated 
by  rounding  errors.  This  weakness  is  a consequence  of  its  use  of  Levinson’s 
recursive  scheme,  which  has  previously  been  reported  to  be  numerically 
unstable  in  some  circumstances  by  Box  and  Jenkins  (1976).  In  contrast, 
our  algorithm  employs  a numerically  stable  technique  to  generate  successive 
systems  of  normal  equations.  The  corresponding  parameters  are  determined 
by  Cholesky's  method,  which  is  known  to  possess  remarkable  numerical 
stability  (e.g.  see  Martin,  Peters,  and  Wilkinson  (1971)). 

Estimation  of  AR  parameters  by  LS  methods  has  previously  been  unpopular 
in  MEM  applications  because  (a)  it  involves  more  computational  effort  than 
Burg's  algorithm,  (b)  the  parameters  may  give  rise  to  a filter  which  is 
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not  minimum  phase,  (c)  the  calculated  values  for  the  parameters  may  "blow 


up”,  and  (d)  the  calculated  value  for  the  residual  sum  of  squares  is 
often  inaccurate  (e.q.  it  may  be  negative!). 

Although  our  algorithm  certainly  involves  more  arithmetic  operations 
than  Burg's  algorithm,  it  is  an  efficient  LS  method  for  estimating 
reasonable  numbers  of  parametets;  later  we  refer  to  some  alternative  LS 
methods  for  problems  involving  larger  numbers  of  patametets.  Secondly, 
if  the  minimum  phase  condition  is  required  in  a given  appll  it  ion  (for 
example,  see  Claerbout  (197b)),  a sensible  proceduie  would  be  to  teflect 
any  poles  resulting  from  the  LS  method  back  inside  the  unit  ciicle.  (Atal 
and  Hanauer  (1971)  discuss  such  a technique,  which  preserves  the  amplitude 
response  of  the  filter.)  Finally,  the  problems  reported  in  (c)  and  (d) 
above  have  been  virtually  eliminated  in  our  algorithm,  through  careful 
attention  being  paid  to  the  computations  involved. 


2.  Statement  of  the  problems 


Given  a time  series  x.,x „,...,x  , let  us  assume  that  x can  be 

1 2 n t 

A 

estimated  by  x^  , where 


m 

V r.  A 

< “ > a,  x , 

t i t-i 

j“l  J J 


for  t « m+1 ,m+2 , . . . ,n  . 


For  given  values  of  the  parameters  (or  filter  coefficients)  a^  , the 
nonrecursive  digital  filter  (1)  provides  a forward  prediction  of  x^  . 

A A 

Putting  et  • xfc  - xfc  and  adopting  the  LS  criterion,  the  parameters 

A r*  A 2 

a are  determined  by  minimizing  the  residual  sum  of  squares  i e . 

1 t-mU 


(Kquivulently , a US  solution  is  sought  for  the  AR  scheme  of  order  m , 
defined  by 
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(2) 


IM 

r»  A A 

•'  ‘ ,1  *i  **->  * ** ' 


Conversely,  suppose  that 


x - £ a.  x 

1 )-l  1 '*3 


for  t » m+ 1 , m+2 , . . . , n .) 


can  be  estimated  by  x^  , where 


for  t - l,2,...,n-m  . 


The  filter  (2)  provides  a backward  prediction  of  x^  , and  its  parameters 


n-m 


— r *2  “ "" 

can  be  determined  by  minimizing  ) e^  , where  - x^  - x^  . 

j t t t t 

It  is  convenient  to  examine  these  problems  in  the  context  of  solving 
systems  of  linear  algebraic  equations.  Thus,  the  forward  prediction 
problem  (1)  is  described  by  the  (n-m)  * m matrix  equation 


(3) 
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X 

n 

If  n > 2m  this  system  is  over deter mined  and  no  a generally  exists  which 
satisfies  (3),  so  it  is  useful  to  define  a residual  vector  e ■ y - ft  a . 

Me  m»  M* 

A LS  solution  a#  to  (3)  is  then  defined  as  any  a which  minimizes  the 

A ATA 

residual  sum  of  squares  S - S(a)  ■ e e , where  the  superscript  T denotes 
the  transpose  operation. 


It  is  well  known  that  a LS  solution  exists  for  any  overdete mined 


It 


system  of  equations,  and  it  is  unique  if  the  matrix  for  the  system  is 
of  full  rank  (■  m , for  problem  (3),  assuming  that  n > 2m).  The  minimum 


value  of  the  residual  sum  of  squares  S is  achieved  when  all  its  first 


partial  derivatives  are  zero  (i.e.  3S/3a^  ■ 0 in  the  case  of  (3)) 


Applying  this  condition  to  (3),  it  follows  that  S is  minimized  when 


AT*  a 

X e ■ 0 , and  so  a#  must  satisfy  the  equation 


(4) 


AT  A A A 

X <y  - a a>  - 2 


The  LS  solution  to  (3)  is  therefore  characterized  by  the  m * m system 

t 


of  normal  equations 


(5) 


ftTA  a ^ta 


X X a. 


For  the  remainder  of  this  paper  we  adopt  a more  concise  notation  and  refer 


to  a4  as  the  solution  to  an  m * m system  of  equations 
(6) 


A A 

5 5. 


A 

S * 


A A , ATA  A , A A A T ATA 

where  R ■ [r  .1  - X X and  s - (s, , s . . . . , s 1 ■*  X y . 

~ 1,3  ~~  - 12  m ~ * 


Similarly,  the  backward  prediction  problem  (2)  is  described  by  the 
(n-m)  * m matrix  equation 


(7) 


X a 


X2 

X3 

X3  Xm+ 1 " 

X4  Xm*2 

"V 

*2 

"X1 

X2 

where  X ■ 

X , 

n-m+ 1 

X ....  X 

n-m*2  n 

# « - 

a 

m 
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l 

x • • • 

3 

1 

3 

1 

+ A A 

From  (4)  we  see  that  the  solution  a„  yields  a residual  vector  e4 


that  is  normal  to  the  column  space  of  X ; hence  the  terminology. 


b 


The  LS  solution  a#  to  (7)  is  characterized  by  the  m * m system  of 
normal  equations 


(tt) 


R «#  “ s , 


where  R - |r.  “XX  and  s ■ Is ^ , s^ • • • • « s^l  “ X ^ . 

Now  in  Burg's  algorithm  for  MEM  it  is  assumed  that  can  be 

estimated  by  a weiqhted  sum  of  m previous  observations  and  a weighted 
sum  of  m future  observations,  using  the  same  weights  a^  in  both 
directions.  Hence  the  MEM  problem  is  described  by  the  2(n-m)  * m 
matrix  equation 


- ,T  -T- 


(9) 


X a 


where  X =* 


A 

X 


, and  y 


A 

y 


The  LS  solution  a#  to  (9)  satisfies  the  m * m system  of  normal 
equations 


(10) 


where  R - lr 


5 S.  3 2 < 


,]  * X X and  s » [s,s,...,s  ] - X y . (If  2n  > 3m 

• j ~ 12  m ~~ 


the  system  (9)  is  overdetermined.  Notice  also  that  we  distinguish 'between 

the  forward,  backward,  and  forward  and  backward  prediction  problems  by 
»A« 

means  of  a superscript,  superscript,  or  no  superscript, 

respectively. ) 
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In  view  of  the  structure  of  X and  ^ , (10)  can  be  rewritten  as 
(11)  (R  + R)a#  * s ♦ s , 

and  comparing  (6) , (8)  and  (11)  it  follows  that  the  MEM  parameter  vector 
aA  is  not  related  in  any  simple  fashion  to  a#  and  a#  . Thus,  the 
normal  equations  (6) , (8)  and  (10)  pertain  to  three  distinct  AR  schemes 
of  order  m , and  our  FORTRAN  program  can  be  directed  to  compute  any  one 

A 

of  the  corresponding  LS  solutions  a#  , or  a#  . These  are  the  three 
variants  of  the  linear  prediction  problem  referred  to  in  the  first 
paragraph  of  §1. 

2-  Thu  algorithm 

Let  us  begin  this  section  by  summarizing  very  briefly  the  current 

state-of-the-art  concerning  algorithms  for  linear  LS  computations.  (The 

book  by  Lawson  and  Hanson  (1974)  provides  a very  detailed  coverage  of  this 

topic,  complete  with  FORTRAN  programs) . 

For  an  overdetermined  system  of  N equations  in  M unknown  para- 

2 t 

meters,  forming  the  normal  equations  requires  NM  /2  operations  and 
solving  them  (by  Cholesky's  method)  requires  M^/6  operations.  However, 
numerical  analysts  generally  prefer  to  avoid  this  approach,  since  forming 
the  normal  equations  considerably  worsens  the  numerical  condition  of  any 
LS  problem  for  which  the  corresponding  over determined  system  of  equations 
is  itself  ill-conditioned. 


+An  operation  is  a multiplication  or  division  plus  an  addition  or  sub- 
traction, and  when  comparing  operation  counts  only  the  highest  powers 
of  N and  M are  given. 


It  is  usually  advocated  instead  that  the  LS  solution  be  computed 


directly  from  the  N x M overdetermined  system  by  an  oithogonalization 

method.  Some  populai  algorithms  of  this  type  are  the  modified  Gram- 

2 

Schmidt  process  (NM~  operations) , Householder  triangular ization 
2 3 

(NM  - M /3  operations) , and  singular  value  decomposition  (at  least 
2 3 

NM  + 5M  operations) . In  a given  precision  floating-point  arithmetic 

these  orthogonalization  methods  successfully  process  a wider  class  of 

LS  problems  than  the  normal  equations  algorithm.  Indeed,  whenever 

single  precision  arithmetic  is  just  adequate  for  computing  the  parameters 

in  a particular  LS  problem  by  (say)  Householder's  method,  the  normal 

( 

equations  must  typically  be  fyrmed  and  solved  in  double  precision  arith- 
metic to  calculate  these  parameters  to  comparable  accuracy. 

Nevertheless,  when  N >>  M (which  is  often  the  case  in  practice)  the 
normal  equations  method  involves  only  about  half  as  many  operations  as  an 
orthogonalization  method.  In  addition,  if  the  normal  equations  can  be 
formed  directly  from  the  data  rather  than  from  an  overdetermined  system, 
then  the  normal  equations  algorithm  requires  only  about  M(M+l)/2  storage 
locations.  Since  orthogot^&lization  methods  require  more  than  NM  storage 
locations  (although  this  can  be  reduced  by  more  sophisticated  programming 
and  some  extra  computation) , the  normal  equations  method  can  offer  worth- 
while savings  in  both  computer  storage  and  time  requirements.  (The  savings 
are  less  significant  when  ill-conditioning  dictates  that  the  normal  equations 
be  formed  and  solved  in  a higher  precision  arithmetic  than  is  required  by 
an  orthogonalization  method.) 

For  the  applications  in  which  we  are  primarily  interested  it  is 
frequently  the  case  that  N > 10M  , and  someti >es  N > 100M  . In  the 


A 
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notation  of  §2,  N = n-m  and  M = m for  (3)  and  (7),  whereas  N = 2(n-m) 
and  M = m in  .-he  MEM  model  (9).  Furthermore,  by  taking  advantage  of 
the  special  structure  of  (6),  (8)  or  (10),  the  normal  equations  for  any 
one  of  these  AR  schemes  can  be  obtained  directly  from  the  given  time 
series  in  [n(m+l)  + 0(m3)]  operations.  We  have  therefore  adopted  the 
normal  equations  approach  for  solving  the  LS  problems  of  this  report. 

In  our  discussion  so  far  it  has  been  tacitly  assumed  that  the  number 
of  parameters  m is  known  a priori,  although  in  practice  this  may  not 
be  true.  Indeed,  choosing  an  appropriate  order  for  an  AR  scheme  is  one 
of  the  most  difficult  tasks  in  time  series  modelling.  Our  algorithm  allows 
this  choice  to  be  made  dynamically  in  a computationally  efficient  manner, 
as  is  explained  later  in  this  section. 

Before  stating  the  general  algorithm,  we  first  show  for  the  normal 
equations  (6)  how  R^  and  s^  (corresponding  to  m = 3)  can  be  obtained 
from  R2  and  s2  (corresponding  to  m = 2) . Suppressing  the  elements 
above  the  diagonal  (since  any  normal  equations  matrix  is  symmetric)  we 
have 


' M2) 
1,1 

• 

“ n-2 

Xt+1  Xt+1 

A(2) 
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n-2  n-2 
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t t+1  t t 
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Notice  that  the  leading  2 * 2 submatrix  of  R^  can  be  obtained  from 

A 

R^  by  subtracting  the  first  term  in  each  summation  (or  inner-product)  which 
defines  an  element  of  R^  , i . e . 

■ K2.\  ■ x3-i  “3-j  • 1 s 1 s 1 s 2 • 

A 

(Equivalently,  the  leading  2*2  submatrix  of  R^  can  be  expressed  as 

Rj  ~ j^2  ' where  the  matrix  ft.,  - lx^  x^]7  fX2  xj^  represents  a rank  one 

modification  to  R^  .)  Also,  the  last  row  of  (apart  from  its  first 

element)  can  be  obtained  from  the  last  row  of  6^  as  follows, 

‘ - *n-J  Vl-j  ' for  2 s 1 s 3 • 


and 


M3)  a(2) 

r3,l  " S2  “ Xn-2  Xn 


Finally,  s^  (apart  from  its  last  element)  can  be  obtained  from  as 


follows 


1 


A ( »)  A(2) 

S a U ” X X 

‘ i ‘ i j-i  3 


for  1 £ i £ 2 , 


...  n—  3 

' ( 3)  i • 

>3  “ L \ xt  + 3 • 


Thus  only  one  inner-product  (the  element  sj^  ) has  to  be  calculated 

afresh  when  we  derive  the  normal  equations  for  m **  3 from  the  normal 

equations  for  m ■ 2 } each  of  the  remaining  elements  of  R^  and  s^ 

is  obtained  at  the  cost  of  one  operation. 

This  scheme  generalises  so  that  6 . and  s , can  be  obtained 

~mi  1 ~m+ 1 

from  R and  s , with  only  the  element  s^m*^  requiring  more  than 
~m  "m  m+1 

one  operation.  When  a backward  prediction  is  needed,  a similar  scheme 

is  available  to  generate  R , and  a , from  R and  s . Further- 

mi»*1  -'inti  •'TO  ~m 

more,  the  one  element  s ' which  must  be  evaluated  as  an  inner- 

m*  1 

product  is  equal  to  s ^ ^ ' , and  so  even  when  the  relationships  (1) 

m*  1 

and  (2)  hold  simultaneously  only  one  inner-product  is  calculated  when 
m increases  to  mtl  . 

Let  us  now  state  the  general  algorithm  which  allows  any  one  of  the 
three  (m+1)  * (mtl)  systems  of  normal  equations  to  be  generated  effi- 
ciently from  the  corresponding  m a m system  of  normal  equations.  The 
algorithm  is  arranged  so  as  to  avoid  the  use  of  intermediate  storage, 
and  since  the  matrices  involved  are  symmetric  the  above- di agonal  elements 


are  never  required. 
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Algorithm  F:  Forward  prediction  problem  (6) 


A 

s - x 


m+1,1  m n-m  n 


for  j - 2, . . . ,m+l  do 


A A 

I‘  f — x X 

m+l,j  m,j-l  n-m  n+l-j 


for  i « 1 , . . . ,m  do 


for  j - 1 , . • . , i do 


A A 

r . , «-  r , . - x , . x . 

1,3  m+l-i  m+1- 3 


for  i *»  1 , . . . , m do 


A A 

s , «-  S , - X , . X 

i i nwl-i  m+1 


m+1 


n-m- 1 

+ l \ *m+l+t 
t-1 


Algorithm  B:  Backward  prediction  problem  (8) 


r . , ♦ s - x,  x 
m+1,1  m 1 m+l 


for  j - 2 , . . . ,m+l  do 


r 4-  r — x x 

m+l,j  m,j-l  j m+1 


for  i - 1 , . . . ,m  do 

for  j - 1, . . . , i do 


r.  . r,  . - x .x 

i,j  i,j  n-m+i  n-m+j 


for  i - 1 , . . . ,m  do 


s,  +•  s - x x 
i i n-m  n-m+i 


m+1 


n-m-1 

l 

t-1 


X X 

t m+l+t 


Algor i t hm  FAB : Forward  and  Backward  predict  ion  problem  (10) 


r , . a - x x - x,  x , 
m+1,1  m n-m  n 1 m+1 


for  j “ 2, • . . , m+1  do 


rm+l,j  * 'm.j-l  " Xn-m  Xn+l-j  “ Xj  “mil 


for  i - 1 , . . . ,m  do 


for  j » 1 , . . . , i do 


*i,j  ri,j  Xmtl-i  Xm*l-j  Xn-m+i  Xn-m+j 


for  i - 1 , .... m do 


s . ♦ s . - x , , x , - x x 
i i mtl-i  m+1  n-m  n-m-»i 


n-m- 1 

s + 2 x V x , 

m+1  t L , t m-tlH 


All  three  algorithms  produce  the  normal  equations  of  order  m+1  by 

first  forming  the  (m+l)th  row,  then  overwriting  the  coefficient  matrix 

and  right-hand  side  vector  of  order  m , and  then  calculating  the  required 

inner-product.  For  Algorithm  F or  Algorithm  H this  step  involves 
2 

(m  +5mt2)/2  operations  plus  n-m- l operations  for  the  inner-product, 

2 

while  Algorithm  FAB  requires  m + Sm-*2  operations  plus  n-m- 1 operations 
for  the  inner-product. 

Starting  at  ml  , and  incrementing  the  value  ot  m by  1 at 

A A 

each  step,  we  can  form  R , R , or  R together  with  s , s , or  s 

for  any  desired  value  of  m . When  ml  there  are  two  inner-product 

| 

calculations  requited  , and  so  the  total  number  of  operations  involved 


t A 

e.g.  in  problem  (t>)  we  calculate  r 


;:i  • T v.  - ■ 1 vm 
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in  forming  the  normal  equations  of  all  orders  from  1 through  m is 


m3  m2  2m 

In(mtl)  — + — - — — - 21  for  Algorithm  F or  Algorithm  B,  and 
6 2 3 


m3  3m2  5m 

(n(m+l)  + — + ■— ---  - — - 3]  for  Algorithm  FAB.  (As  a comparison, 

i / n 


3 2 6 

if  the  normal  equations  were  all  derived  independently  with  inner- 


products  calculated  from  first  principles,  the  operation  count  in  forming 

3 


problem  (6)  for  orders  1 through  m would  be  n i m‘  + — 1 

6 6 


4 * -,2 

. m Jm  /m  m.  . ...  . _ * . . 

[ — + — — + — — + — ] ; when  n = 100  and  m = 10  this  approach 
8 4 8 4 


involves  25,410  operations,  whereas  Algorithm  F requires  only  1,308 
operations . ) 

Assuming  that  an  appropriate  value  for  m is  available,  the  relevant 


normal  equations  of  order  m can  be  formed  by  one  of  the  Algorithms  F,  B, 

3 


or  FAB,  and  then  solved  by  Cholesky's  method  (m  /6  operations).  However, 
let  us  suppose  that  the  user  is  not  sure  what  value  for  m is  appropriate. 


but  that  he  can  specify  values  m and  m,  _ such  that  m _ _ i m £ m, 

start  last  start  last 


In  this  event  a sensible  algorithm  would  not  solve  the  normal  equations  until 


the  value  of  m has  been  increased  to  m , even  though  the  normal 

start 


equations  are  formed  for  m=l,2,...,n> 

start 


Our  algorithm  employs 
this  strategy,  and  so  if  the  user  knows  the  "correct"  value  for  m before- 


hand he  simply  sets  m *-  = m,  However,  typical  applications  of 

StdTt  IdSt 


our  algorithm  call  for  a sequence  of  normal  equations  to  be  solved,  and 


hence  beginning  at  m = mstart  we  generate  a duplicate  copy  of  the  current 


normal  equations  since  these  are  subsequently  destroyed  by  the  Cholesky 
factorization. 

The  problem  of  choosing  an  appropriate  value  for  the  number  of  para- 
meters m depends  upon  the  particular  application  involved.  (The  choice 
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might  also  depend  on  the  amount  of  inherent  error  in  the  time  series, 
and  on  the  accuracy  of  the  computing  environment.)  For  MEM  applications 
we  often  use  a criterion  due  to  Akaike  to  determine  the  order  m of  the 
AR  scheme;  see  Appendix  A for  a description  of  this  criterion.  The 
program  given  in  Appendix  E includes  a subroutine  which  implements 
Akaike' s scheme,  and  so  users  who  wish  to  determine  m differently  can 
simply  replace  this  subroutine  with  one  of  their  own  choosing. 

Most  of  our  digital  signal  processing  applications  involve  the  use 
of  minicomputers  which  operate  in  fixed  precision  floating-point  arith- 
metic. Hence,  the  single  precision  program  in  Appendix  E is  designed 
to  avoid  any  dependence  on  higher  precision  program  segments.  Conse- 
quently it  is  imperative  to  take  extra  care  when  performing  any  floating- 
point calculations  where  serious  loss  of  significance  can  occur.  Compu- 
tational experience  with  our  algorithm  confirms  that  significance  losses 
can  be  troublesome  when  calculating  inner-products  (particularly  when 
n is  large),  and  also  in  computing  the  residual  sum  of  squares.  These 
two  computations  are  discussed  in  Appendices  C and  D respectively,  and 
we  provide  accurate  subroutines  for  both  calculations  in  Appendix  E. 

Finally,  the  computational  efficiency  of  our  algorithm  has  been 
optimized  in  its  FORTRAN  implementation  by  avoiding  the  use  of  two- 
dimensional  arrays.  Hence  Appendix  B contains  an  alternative  description 
of  Algorithm  FAB  in  vector  notation. 

4_.  Further  details  and  alternative  algorithms 

The  algorithm  described  in  § 3 generates  each  successive  system  of 
normal  equations  from  the  preceding  system,  but  each  set  of  equations 


i 


lb 


is  solved  independently  (i.e.  without  reference  to  the  cholesky 

solution  of  previous  systems).  Now  Cholesky's  method  requires  about 
J 2 t 

m /b  ♦ tm  /2  t 0(m)  operations  to  solve  any  m * m symmetric 

positive  definite  * system  ot  lineal  equations.  Hence,  if  the  normal 

equations  are  solved  independently  for  a_H  orders  from  1 through  m , 

4 3 2 

Cholesky's  method  uses  a total  of  about  m /24  ♦ 7m  /12  + 0 (m  ) 

operations.  Nevertheless,  in  out  typical  applications  the  work  required 

by  Algorithm  K to  form  all  the  normal  equations  is  still  more  than  the 

2 

work  of  solving  them;  tor  example,  tins  is  the  case  when  n S m and 
ms  15  . 

In  this  section  we  tefet  briefly  to  some  algorithms  which  do  take 

advantage  of  previous  solutions.  Wo  restrict  our  remarks  primarily  to 

the  forward  prediction  problem  (6),  but  analogous  alternative  methods 

exist  tor  the  backward  prediction  problem  (H)  and  for  the  MEM  problem 

(10).  Although  the  method  of  S3  is  usually  more  efficient  than  these 

alternative  schemes  when  ms  15  , they  become  more  attractive  for 

larger  values  ot  m . However,  they  also  require  more  programming 

effort  than  is  involved  m Appendix  E,  particularly  when  appropriate 

techniques  to  pieserve  numei ical  accuracy  are  incorporated. 

Let  us  begin  by  restating  the  relationship  that  exists  between 

1 


A rA(mtl) 

~mt 1 L 1 » i 


n-m-  l 

l * 

t-1 


m* 1- i tt  ’ mt 1- j tt 


and 


a Ta  (m) 

R ■ r . , «* 

-m  l l,  j_ 


t+ 


The  0(m)  term  in  this  operation  count  includes  m square  root  reciprocals, 
although  an  alternative  version  exists  in  which  square  roots  are  avoided. 

Any  normal  equations  matrix  is  positive  definite  provided  that  the  cor- 
responding overdetermined  system  ot  equations  is  of  full  rank. 


i.V 

5$ 


I 


y x x . Letting  A be  an  m x m matrix  with  elements 

m-i+t  m-j+t  ~m 

t»l 

U J 

a (m)  , ^ fA(m+l)  a (m+1) 

6 , and  adopting  the  more  concise  notation  p “ r , r -<•••* 

i,j  r ~m  Lm+1,1  m+1, 2 

Mm+1)  l an<j  „ Mmtl)  we  see  that  Algorithm  F provides  a 

m+l,mj  m+1  m+1, m+1 

constructive  proof  of  the  following  result. 

Theorem  1.  For  the  forward  prediction  problem  (6) 


A A 

R - A 
~m  ~m 


where  Sf15*!  ■ x . x 

i,j  m+l-i  m+l-j 


Theorem  l shows  that  R . (apart  from  its  extra  row  and  column)  can  be 

■m+1 

A A 

obtained  by  subtracting  a rank  one  matrix  A from  R , where 


lx  ,x 
m m- 1 


,x  ] |x  ,x 
1 m m- 1 


When  solving  the  m x m normal  equations  (6)  by  Cholesky's  method, 

the  first  step  (o(m')  operations)  is  to  obtain  a factorization 
a at  a a , , 

L L of  R , where  the  matrix  1.  is  lower  triangular.  The  second  step 
~nv~m  "m  ~m 

I 2 \ a at  a A 

lo(m  ) operations  obtains  the  solution  to  L L a = s , by  solving 

v ' 'm~m  ~m  •'in 

AAA  AT  A A 

Lb  = s and  then  L a b . Thus,  whenever  problem  ( fc> ) is  to  be 

~m'm  ~m  ~ '-m  'm 

solved  for  many  successive  values  of  m it  is  worthwhile  improving  the 
efficiency  of  the  f irst  step. 

In  view  of  Theorem  1,  this  may  be  accomplished  by  obtaining  the 

Cholesky  factorization  L , L , of  R , from  the  previous  factorization 

~m+l  m+1  ~m+l 

L I.1"  of  R . Methods  for  updating  the  factors  of  a symmetric  positive 
~m~m  ~m 

definite  matrix  following  a rank  one  modification  are  described  by  Gill, 


Golub,  Murray  and  Saunders  (1974),  who  also  explain  how  these  methods  can 


is 


be  used  to  process  the  extra  row  and  column  (exhibited  in  Theorem  1)  of 

A 

R . (This  special  case  relies  on  the  fact  that  the  rectangular  matrix 

A A 

X can  be  obtained  tiom  X by  deleting  one  of  its  rows  and  adding  a 

~m+ 1 'm 

new  column.)  When  modifying  Cholesky  factorizations  the  problem  of 
maintaining  positive  definiteness,  in  the  presence  of  rounding  error, 
must  be  handled  carefully  (see  also  Fletcher  and  Powell  (1974)),  but 

numerically  stable  methods  are  available.  The  most  efficient  of  these 

2 A 
require  cm  * 0(m)  operations  to  obtain  the  factors  of  R from 

A 

those  of  R , where  3 < c < 5 . (In  contrast,  it  only  requires 

'TO 

2 A 

m + 0(m)  operations  to  obtain  the  factors  of  R from  those  of 

'TO 

A 

R . In  our  applications  however,  it  is  more  convenient  to  increase 
"Tn+ 1 

m than  to  decrease  it.) 

For  the  MEM  problem  (10)  the  matrix  R , (apart  from  its  extra  row 

~m+l 

and  column)  can  be  obtained  by  subtracting  a rank  two  matrix  A from 

~m 

A 

R . The  matrix  A is  the  sum  of  A and  its  counterpart  A from 

~m  ~m  ~m  ~m 

Algorithm  B,  which  has  the  form 


A = [x  ,x  ,...,x  1 lx  .,x  ,x  ] 

~m  n-m+1  n-m+2  n n-m+1  n-m+2  n 


The  factors  of  R , can  be  obtained  from  the  factors  of  R either  by 
~m+ 1 '-m 

applying  one  of  the  above  rank  one  methods  twice,  or  by  using  a rank  two 

method  (see  Bennett  (1965)  for  example) . 

Thus,  solving  the  normal  equations  (6) , (8) , or  (10)  for  all  orders 

from  1 through  m by  methods  involving  modifications  to  Cholesky 

3 2 

factorizations  requires  (c'+l)m  /3  + 0 (m  ) operations,  where  c'  = c 
for  (6)  or  (8)  and  c'  = 2c  for  (10) . For  moderate  values  of  m these 
methods  therefore  offer  no  gain  in  efficiency  for  the  algorithm  of  §3. 
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However,  in  a very  recent  paper  by  Morf,  Dickinson,  Kailath  and 

2 

Vieira  (1977)  this  operation  count  has  been  reduced  to  7m  + 15m  - 7 
for  problems  (6)  or  (8) : they  do  not  consider  the  MEM  problem  (10) . 

This  order-of-magnitude  improvement  in  the  effort  required  to  solve  all 
m normal  equations  is  made  possible  because  the  matrix 

R = (y  | X]  ly  | X] 

is  "close"  to  Toeplitz  form,  in  a sense  made  precise  by  Friedlander, 

Morf,  Kailath  and  Ljung  (1977).  Their  algorithm  also  takes  advantage 

of  the  product-form  structure  of  & , since  it  is  the  fact  that  [y  | X] 

is  a rectangular  Toeplitz  matrix  which  permits  the  equations  to  be  solved 

so  efficiently  by  an  extension  of  Levinson's  recursive  scheme. 

In  view  of  the  fact  that  Levinson's  original  algorithm  for  (square) 

Toeplitz  matrices  is  numerically  unstable  in  some  circumstances,  it  seems 

likely  that  this  extended  algorithm  also  suffers  from  occasional  rounding 

error  difficulties.  Indeed,  because  of  their  high  speed  and  low  storage 

requirements,  it  would  clearly  be  advantageous  to  develop  numerically 

stable  versions  of  both  the  original  and  extended  Levinson  algorithms. 

, 2 

(Hopefully,  this  could  be  accomplished  while  retaining  the  0(m  ) 
operation  count.) 

Finally,  although  the  above  summary  of  alternative  algorithms  gives 
some  prominence  to  the  number  of  operations  involved  in  solving  all  m 
normal  equations,  the  additional  n(m+l)  operations  required  by  all 
the  methods  often  dominates  the  total  operation  count.  Thus,  when 
n >>  m , the  gains  in  effioiency  offered  by  these  alternative  methods 
may  be  more  apparent  than  real.( 


20 


5.  The  max  imam  ent  ropy  spectrum 

Our  primary  motivation  for  developing  the  program  in  Appendix  E 
was  to  determine  the  AR  parameters  tor  MEM  spectral  analysis.  Hence, 
for  the  sake  of  completeness,  in  this  section  we  provide  a brief 
heuristic  description  of  the  spectrum  calculation,  followed  by  a 
numerical  example  which  compares  spectral  estimates  obtained  by  our 
algorithm  and  by  Burg's  method.  (The  theory  of  MEM  and  its  relation- 
ship to  AR  processes  has  been  described  in  numerous  other  reports;  in 
particular,  see  Ulrych  and  Bishop  (1975).) 

Given  a time  series  x,,x  . ...,x  , for  which  At  is  the  uniform 

12  n 

sampling  interval,  the  first  step  is  to  apply  Algorithm  FAB  to  obtain 


the  LS  solution  a.  = la  .a...., a ] and  its  associated  minimum 

12  m 


residual  sum  of  squares  S = S (a  ) . This  yields  a recursive 

m m 


predictive  filter  for  the  signal,  whose  mean  square  residual  P is 

m 


given  by 


(12) 


p = S 

m 2(n-m)  m 


If  the  fit  to  the  time  series  is  good,  the  residuals  e^  should  have 


small  correlation,  resalting  essentially  in  a white  noise  process. 


When  the  MEM  parameters  are  used  for  forward  prediction, 


the  residual  e is  defined  as 
t 


In  the  case  of  Algorithm  F (or  Algorithm  B)  the  mean  square  residual 
is  obtained  by  dividing  the  residual  sum  of  squares  by  (n-m)  . 


I 


J 


The  z transform  of  the  difference  equation  (13)  yields 


(14)  E(z)/X(z)  = 1 - 7 a*  z”3 

j-1  3 

where  z ^ is  the  unit  delay  operator.  When  the  time  series  of 

interest  is  input  to  the  filter  defined  by  (13),  it  produces  (approximately) 

a white  noise  output.  Thus,  if  a white  noise  sequence  is  applied  to 

the  inverse  of  this  filter  the  resulting  time  series  has  a spectrum 

similar  to  that  of  the  original.  When  the  residual  is  white,  E(z) 

is  a constant  satisfying 

(15)  |E(z) |2  = P At 

' m 

where  P is  defined  by  (12).  Solving  (14)  for  X(z)  , and  taking 
m 

the  square  of  the  modulus,  we  obtain  the  expression 


(16) 


P At 

m 

V * -•tiufjAt 
> a . e 

j.l  3 


for  the  spectrum  P(f)  at  frequency  f . (The  Fourier  transform  is 

obtained  by  evaluating  the  z transform  on  the  unit  circle,  by  sub- 

...  . < 2 ti  f At 

stituting  z = e . ) 

All  that  remains  is  to  evaluate  (16)  for  as  many  values  of  f as 
desired,  between  0 and  the  Nyquist  (folding)  frequency  l/(2At)  . 


When  periodic  signals  are  present  the  spectrum  must  be  evaluated  at 


22 

close ly— spaced  frequencies  in  order  to  obtain  a satisfactory  representation 
of  its  peaks. 

We  conclude  this  section  by  illustrating  the  superiority  of  LS 
solutions  to  the  MEM  problem  for  short  sample  sinusoidal  data.  Our  test 
signal  was  formed  by  summing  a 0.0J  Hz  and  a 0.2  Hz  sine  wave,  generated 
in  single  precision  arithmetic  (approximately  6 significant  digits)  and 
sampled  10  times  per  second.  Subroutine  FABNE  (from  Appendix  E)  and 
Burg's  method  (as  described  by  Andersen  (1974))  were  used  to  obtain  MEM 
spectral  estimates  by  processing  n = 75  data  points.  The  resulting  spectra 
are  shown  in  Figure  1. 

The  number  of  cycles  available  for  analysis  are  1.5  and  0.23  for  the 
0.2  Hz  and  0.03  Hz  signals,  respectively.  The  spectra  are  shown  one 

above  the  other  for  increasing  values  of  m (the  number  of  coefficients). 

4 

The  vertical  scale  is  80  dB  per  division  (a  factor  of  10  in  power) , 
and  successive  spectra  are  shifted  up  by  80  dB  . We  have  drawn  in  the 
test  signal  frequencies  for  convenience. 

Since  the  purpose  of  this  example  is  to  compare  LS  solutions  with 
those  obtained  from  Burg's  method,  we  used  double  precision  versions  of 
both  subroutine  FABNE  and  Burg's  algorithm  in  order  to  minimize  rounding 
error  effects^.  Notice  that  the  Burg  algorithm  causes  frequency  splitting, 
and  that  it  does  not  converge  as  m increases.  In  contrast,  subroutine 


In  single  precision  arithmetic  subroutine  FABNE  produces  correct  results 
out  to  m = 6 , but  at  m = 7 it  returns  an  error  code  indicating  that 
positive  definiteness  has  been  lost. 


Comparison  of  BURG  and  FABNE  spectra  for  a 0.2  and  0.03  Hz  sine  wave 
Spectra  are  plotted  for  increasing  numbers  of  coefficients. 
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PAHNK  converges  to  the  correct  answer  for  the  0.2  Hi:  sine  wave  when 
m » 3 , and  the  0.03  Hz  signal  is  clearly  evident  when  m * 16  . 

Finally,  in  fairness  to  Burg's  method,  our  experience  suggests  that 
when  many  cycles  of  a sinusoid  are  available  for  analysis  (say,  15  or 
more  cycles)  the  spectral  estimates  obtained  from  both  methods  are 
usually  quite  similar. 


Acknowledgements 

We  are  indebted  to  Dr.  J.E.  Lokken  for  his  initial  support  and 
continued  interest  in  out  work,  to  Mr.  K.U.  Wilson  for  programming 
assistance,  and  to  the  Chief  of  Research  and  Development,  Department 
of  National  Defence,  tor  financial  assistance.  We  have  also  benefitted 
greatly  from  discussions  (on  matters  related  directly  and  indirectly 
to  the  contents  of  this  report)  with  Professors  T.J.  Ulrych,  J.F.  Claerbout, 
M.  Morf,  T.  Kailath,  G.H.  Golub,  L.M.  Delves,  W.K.  Hastings,  R.R.  Davidson, 


and  M.J.D.  Powell. 


References 


Akaike,  H.,  1969.  "Fitting  autoregressive  models  for  prediction".  Ann. 

Inst.  Stat.  Math.,  vol.  21,  pp.  243-247. 

Andersen,  N.,  1974.  "On  the  calculation  of  filter  coefficients  for 

maximum  entropy  spectral  analysis".  Geophysics,  vol.  39,  pp.  69-72. 
Atal , B.S.  and  Manauer,  S.L.,  1971.  "Speech  analysis  and  synthesis  by 

linear  prediction  of  the  speech  wave".  J.  Acoust . Soc.  Amer.,  vol. 

60,  pp.  637-656. 

Barber,  F.G.  and  Taylor,  J.,  1977.  "A  note  on  free  oscillations  of 

Chedabucto  Bay".  Dept.  Environment.  Mar.  Sci.  Directorate,  Man. 

Rep.  Series  No.  47,  Ottawa,  Ont. 

Bennett,  J.M.,  1965.  "Triangular  factors  of  modified  matrices".  Numer. 
Math.,  vol.  7,  pp.  217-221. 

Box,  G.E.P.  and  Jenkins,  G.M.,  1976.  Time  Series  Analysis.  Holden-Day, 

San  Francisco,  Calif. 

Burg,  J.P.,  1967.  "Maximum  entropy  spectral  analysis".  37th  Ann.  Int. 

Meet.,  Soc.  Explor.  Geophys . , Oklahoma  City,  Ok  la. 

Burg,  J.P.,  1968.  "A  new  analysis  technique  for  time  series  data". 

NATO  Adv.  Study  Inst,  on  Signal  Processing,  Enschede,  Netherlands. 

Burg,  J.P.,  1975.  "Maximum  entropy  spectral  analysis".  Ph.D.  thesis, 
Stanford  University,  Palo  Alto,  calif. 

Chen,  W.Y.  and  Stegen,  G.R.,  1974.  "Experiments  with  maximum  entropy 

power  spectra  of  sinusoids".  J.  Geophys.  Res.,  vol.  79,  pp.  3019-3022. 
Claerbout,  J.F.,  1976.  Fundamentals  of  Geophysical  Data  Processing. 
McGraw-Hill,  New  York. 

Dahlquist,  G.  and  BjBrk,  A.,  1974.  Numerical  Methods.  Prentice-Hall, 


Englewood  Cliffs,  N.J. 


Erickson,  H.E.,  1976.  "Some  experiments  with  maximum  entropy  spectral 
analysis".  23rd  Pacif.  N.W.  Req . Meet.,  Amer.  Geophys.  Union, 

University  ot  Victoria,  Victoria,  B.C. 

T 

Fletcher,  R.  and  Powell,  M.J.D.,  197*1.  "On  the  modification  of  LDL 
factorizations'.  Math.  Comp.,  vol.  28,  pp.  1067-1087. 

Friedlander,  B.,  Morf,  M.,  Kailath,  T.  and  Ljunq,  L. , 1977.  "New 

inversion  formulas  for  matrices  classified  in  terms  of  their  distance 
from  Toeplitz  matrices".  Submitted  for  publication. 

Gill,  P.E.,  Golub,  G.H.,  Murray,  W.  and  Saunders,  M.A.,  1974.  "Methods 

for  modifying  matrix  factorizations" . Math.  Comp.,  vol.  28,  pp.  505-535. 

Jones,  R.H.,  1976.  "Autoregression  order  selection".  Geophysics,  vol.  41, 
pp.  771-773. 

Kahan,  W.,  1972.  "A  survey  of  error  analysis",  pp.  1214-1239,  from 
Information  Processing  71.  Nor th-Hol land,  Amsterdam. 

Lawson,  C.L.  and  Hanson,  R.J.,  1974.  Solving  Least  Squares  Problems. 
Prentice-Hall,  Enqlewood  Cliffs,  N.J. 

Levinson,  N.,  1947.  "The  Wiener  RMS  error  criterion  in  filter  design  and 
prediction".  J.  Math.  Phys.,  vol.  25,  pp.  261-278. 

Martin,  R.S.,  Peters,  G.  and  Wilkinson,  J.H.,  1971.  "Symmetric  decomposition 
of  a positive  definite  matrix",  pp.  9-30,  from  Handbook  for  Automatic 
Computation , vol.  II,  by  Wilkinson,  J.H.  and  Reinsch,  C.  Springer- 
Verlag,  New  York. 

Morf,  M. , Dickinson,  B.,  Kailath,  T.  and  Vieira,  T. , 1977.  "Efficient 

solution  of  covariance  equations  for  linear  prediction".  IEEE  Trans. 
Acoust.,  Speech,  Signal  Processing,  vol.  ASSP-25,  pp.  429-433. 


i 


1 


Ulrych,  T.J.  and  Bishop,  T.N.,  1975.  "Maximum  entropy  spectral  analysis 

and  autoregressive  decomposition".  Rev.  Geophys.,  vol.  13,  pp.  183-200. 
Ulrych,  T.J.  and  Clayton,  R.W.,  1976.  "Time  series  modelling  and  maximum 
entropy".  Phys . Earth  Planet.  Inter.,  vol.  12,  pp.  188-200. 


Appendix  A:  Aka ike 1 s final  ptediction  error (FPE)  criterion 

The  deteinunut ion  of  the  or  do  t of  the  AK  process  is  of  central 
importance  in  MEM  spectral  analysis.  Akaike's  (19b9)  final  prediction 
error  (FPE)  for  an  n- point  time  series  is  given  by 


(Al) 


F<PF 

m 


n ♦ m » 1 p 

n - m - 1 m 
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n 


:J 


where  P denotes  the  mean  square  lesidual  (12)  lor  the  1.S  solution 
m 

a#  to  the  m parameter  problem  (10).  The  "optimal"  AK  order  m is 

then  determined  by  the  minimum  value  of  FPE  for  m « 1,2,...  . 

m 

If  the  mean  is  not  subtracted  fiom  the  data  prior  to  the  calculation  of 
the  AR  parameters,  expression  l A l ) must  be  replaced  by 


( A2 ) 


FPE 

m 


n * m 
n - m m 


Various  opinions  on  the  effectiveness  of  the  FPE  criterion  have 

been  expressed  (e.q.  see  Ulrych  and  Bishop  (1975),  Ulrych  and  Clayton 

(197b),  Jones  (197b),  Barber  and  Taylor  (1977),  and  references  therein), 

but  most  of  this  experience  relates  to  its  apt  lication  to  Burg’s  algorithm. 

Our  experier.ee  with  the  LS  algorithm  of  this  paper  suggests  that  while 

the  FPE  criterion  provides  a useful  automatic  method  of  selecting  an 

appropriate  value  for  m , it  should  be  used  with  caution.  (If  m is 

too  small  the  resulting  spectrum  will  be  too  smooth,  and  if  m is  too 

large  extraneous  peaks  will  occur.) 

For  any  of  the  three  linear  prediction  problems  (6),  (8),  or  (10) 

ths  minimum  residual  sum  of  squares  s decreases  as  m increases. 

m 

Therefore,  even  when  the  statistical  assumptions  underlying  the  FPE  criterion 


do  not  hold,  it  seems  reasonable  to  allow  m to  be  determined  as  the 


smallest  order  for  which  S £ S , . When  n >>  m this  can  be 

m m+1 

accomplished  by  choosing  m to  be  the  first  local  minimum  of  FPE 

m 

Thus,  subroutine  AKSTOP  of  Appendix  E determines  the  "optimal"  m 

for  problems  (6),  (8),  or  (10)  as  the  smallest  integer  in  the  interval 

[m  ,m,  1 such  that  FPE  , > FFe 

start  last  m+1  m 


Appendix  B:  Vector  notation  employed  in  the  FORTRAN  program 

Since  some  FORTRAN  compilers  do  not  process  two-dimensional 
arrays  very  efficiently,  our  program  uses  only  vectors  (i.e.  linear 
arrays)  when  forming  and  solving  the  normal  equations.  Thus,  the 

A 

lower  triangle  of  the  matrix  R (or  R or  R)  is  stored  by  rows 
as  a vector  v with  m(m+l)/2  elements.  For  example,  when  m=4 
we  have : 


r r 

3,2  r3,3 


V2  V3 


V4  V5  V6 


r4,2  r4,3  r4,4j 


Lv  V V V 
7 8 9 10 


For  the  convenience  of  readers  who  require  a statement  of  our 
algorithm  that  is  closely  related  to  the  code  in  Appendix  E,  we  pro- 
vide the  following  description  of  Algorithm  FAB  in  this  vector  notation: 


k + m(m+l)/2  + 1 


V,  +-s-x  X - X.  X , 
k m n-m  n 1 m+1 


j - 1 
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for  i - k+l,...,k+m  do 

j ♦ j ♦ 1 

i l-m- 1 n-m  n+l-j  j m+1 

k ♦ 0 

for  i - 1 , . . . , m do 

for  j = 1 , . . . , i do 


k + k+1 

v.+-v,  - x.,  -x.,.-x  ..x 

k k m+l-i  m+l-j  n-m+i  n-m+3 

for  i « 1 m do 

S . S , “X.,  .x.,  “X  X 

x 1 m+l-i  m+1  n-m  n-m+i 


m+1 


•+  2 * 


n-m- 1 

l 

t=l 


t m+l+t 


Appendix  C:  Accurate  calculation  ot  inner-products 

Inner-product  calculations  occur  frequently  in  scientific  program- 
ming. Indeed,  in  any  system  of  normal  equations  every  element  of  the 
coefficient  matrix  and  the  right-hand  side  vector  is  defined  as  an 
inner-product.  In  a typical  program  an  algorithmic  statement  such  as 

N 

(Cl)  sum  + £ x w 

t-1  t t 

is  usually  implemented  in  floating-point  arithmetic  by  computing 
N 

£ x w from  single  precision  data  x and  w by  forming  double 
t-1  * 

precision  products  and  double  precision  sums,  and  then  assigning  the 
final  result  to  the  single  precision  variable  sum.  However,  for  many 


31 
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i < 


of  our  applications  it  is  necessary  to  avoid  any  dependence  on  higher 
precision  arithmetic,  and  so  the  above  scheme  cannot  be  employed. 

A single  precision  FORTRAN  program  segment  of  the  form 

(C2)  SUM  =0.0 

DO  10  I = 1,  N 
10  SUM  = SUM  + X(I) *W(I) 

can  produce  a calculated  value  for  SUM  which  differs  from  its  true 

N 

value  by  almost  1.06u  * £ |(N+2-t)x  w | , where  u represents  the 

t=l 

roundoff  unit  of  the  floating-point  system  used  (e.g.  see  Dahlquist  and 
Bjttrk  (1974)).  Even  though  the  actual  error  incurred  is  often  consider- 
ably less  than  this  upper  bound,  for  large  values  of  N the  program 

(C2)  does  not  usually  compute  inner-products  accurately  enough  when 

-6 

forming  the  normal  equations.  (For  example,  if  n = 500  , u s 10  , 

|xt(  SI  , and  | w ^ | v 1 , typical  results  from  (C2)  contain  only  about 

four  significant  decimal  figures.) 

We  have  therefore  adopted  the  following  single  precision  program 

N 

segment  when  computing  ) x w for  use  in  forming  the  normal  equations 

t=l  t t 

(6),  (8),  or  (10).  (It  is  based  on  a program  due  to  Kahan  (1972)  for 
N 

computing  £ x , which  has  a very  small  error  bound  that  is  essentially 
t=l  t 

independent  of  N . ) 


SUM 

= 

0. 

0 

C = 

0. 

0 

DO  , 

20 

1 

= 

1. 

N 

Y - 

C 

+ 

X 

(I) 

*W(I) 

T - 

SUM 

+ 

Y 

C - 

(SUM 

l-T) 

♦ Y 

20  SUM 

- 

T 

SUM 

= 

SUM 

•f 

C 

1 


1 


<C3) 
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The  role  of  the  variable  C in  (C3)  is  to  recover  the  information  that 
is  lost  when  adding  the  products  xtwt  * t>ut  notice  that  no  attempt  is 
made  to  increase  the  precision  of  these  products.  Thus  (C3)  simulates 
double  precision  accumulation  of  single  precision  products,  which  is 
usually  satisfactory  when  SUM  is  of  large  magnitude  compared  to  these 
products . 

This  is  almost  always  the  situation  that  occurs  in  our  applications 

where  most  of  the  N products  xtwt  are  positive.  The  reason  for  this 

is  because  the  term  w in  (Cl)  is  replaced  by  x when  forming 

t m+l+t 

the  normal  equations,  and  our  typical  time  series  xj'X2*”'*'xn  have 
few  sign  changes  relative  to  n (i.e.  most  pairs  x ,x  are  of  the 

same  sign  when  n >>  m ) . 

Subroutine  ACCSUM  of  Appendix  E is  based  on  the  program  segment 
(C3) . It  is  used  to  compute  each  inner-product  that  is  calculated  from 
first  principles;  thus  it  is  called  twice  when  m=l  and  once  for  each 
value  of  m = 2,3,...  . Provided  that  SUM  is  large  enough  to  ignore 
the  rounding  errors  made  in  forming  the  products  xtwt  * subroutine 
ACCSUM  produces  a calculated  value  of  SUM  for  which  the  relative 
error  is  bounded  by  5u  . This  bound  applies  to  any  practical  value  of 
N , and  in  our  experience  the  actual  magnitude  of  the  relative  error 
rarely  exceeds  u when  N < 1,000  . 


I 


Appendix  D:  Calculating  the  residual  sum  of  squares 

Given  an  overdetermined  N * M system  of  equations  X a = y , 

T 

where  X =*  [x  .]  and  y = [y  ,y  , — ,y  ] , the  residual  sum  of 

^ m]  ~ * * N 
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squares  S ■ S (a)  is  defined  by 


S (a)  = eTe  «=  (y  - X a)T(y  - X a) 


(Dl) 


N M _ 

I (y.  - I a.x.  .) 

iil  1 j=l  3 x'j 


For  any  given  a , evaluating  S directly  from  the  second  line  of 
(Dl)  involves  NM  operations.  Although  X and  y take  special  forms 
in  the  linear  prediction  problems  (3),  (7),  and  (9),  their  structure 

yields  no  reduction  in  this  operation  count.  (Recall  that  N = n-m  for 
(3)  and  (7),  N = 2 (n-m)  for  (9),  and  M = m in  all  three  problems.) 
Thus,  when  n >>  m , calculating  the  residual  sum  of  squares  via  (Dl) 
for  any  value  of  m requires  almost  as  much  computation  as  forming  the 
normal  equations  of  all  orders  from  1 through  m by  Algorithms  F,  B, 
or  FAB.  Furthermore,  when  solving  a sequence  of  normal  equations  (6), 

(8) , or  (10)  for  successive  values  of  m , S must  be  evaluated  by  (Dl) 
ab  initio  for  each  value  of  m . 

In  view  of  this  it  is  quite  understandable  that  the  following  fast 
method  of  calculating  S is  popular,  particularly  among  engineers  and 
statisticians.  Expanding  the  first  line  of  (Dl)  we  obtain 

T T T T T T 

(D2)  S (a)  = y y - a (X  y)  + a (XXa-Xy)  , 

and  since  the  LS  solution  a#  satisfies  the  normal  equations 

T T 

(D3)  X X a,  » X y , 

substituting  l D 3 ) into  (D2)  gives 


14 


(D4) 


T T T 

SM<'*)  * t X “ V 


Calculating  the  minimum  residual  sum  ot  square,  by  formula  (D4)  involves 

T 

only  N ♦ M operations,  since  the  term  X y was  formed  previously  as 

the  right-hand  side  vector  of  the  normal  equations  (D3) . 

A further  reduction  can  be  made  in  this  operation  count  when  solving 

a sequence  of  normal  equations  (b) , (8),  or  (10)  for  successive  values  of 

T 

m . Specifically,  when  m - 1 the  term  y y in  (D4)  can  be  obtained 

T 

in  two  operations  from  X X . For  example,  in  the  case  of  (6),  when 

ata  ata  2 2 

m = 1 we  calculate  y y as  (X  X - x + x ) . Thereafter,  for  each 

~ ~ ~ 1 n 

ata 

m = 2,3,...  , we  obtain  the  current  value  of  y ^ in  one  operation  via 

ata  ata  2 

the  assignment  statement  y y *•  y y - x . Problem  (8)  is  handled  in  an 

~ ~ ~ ~ m 

analogous  fashion,  while  problem  (10)  requires  two  operations  to  update 
T 

the  value  of  y y for  each  m . Thus,  for  any  one  of  the  three  linear 
prediction  problems,  calculating  the  minimum  residual  sum  of  squares  by 
formula  (D4)  involves  no  more  than  m+2  operations  for  each  successive 
value  of  m . 

Our  program  in  Appendix  E can  be  directed  to  obtain  the  minimum 
residual  sum  of  squares  in  this  efficient  manner,  merely  by  setting  the 
input  parameter  RSCODE  of  subroutine  FADNE  equal  to  0 . Alternatively, 
setting  RSCODE  = 1 causes  this  minimum  to  be  calculated  via  subroutine 
ACCSSQ  , which  employs  formula  (Dl)  rather  than  (D4) . 

As  is  clear  from  the  above  operation  counts,  the  use  of  subroutine 
ACCSSQ  is  expensive  in  terms  of  computer  time,  but  our  computational 
experience  leaves  no  doubt  as  to  the  necessity  for  this  optional  subroutine 
in  practice.  Formula  (D4)  involves  a subtraction  of  two  calculated  terms 
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which  are  almost  always  of  the  same  magnitude,  and  which  frequently  agree 

to  more  than  two  or  three  significant  decimal  figures.  Indeed,  this  has 

to  be  the  case  when  the  minimum  residual  sum  of  squares  is  small  relative 
T 

to  y y , and  a single  precision  implementation  of  formula  (D4)  is  then 
bound  to  suffer  from  significance  losses.  (In  our  opinion,  many  of  the 
reported  difficulties  in  implementing  "optimal”  selection  rules  (such  as 
Akaike's  criterion)  can  be  attributed,  at  least  in  part,  to  careless 
Usage  of  formula  (D4) .) 

Since  our  algorithm  forms  and  solves  the  normal  equations  in  a numeri- 
cally stable  manner,  it  is  not  possible  to  significantly  increase  the 
accuracy  obtainable  from  (D4)  without  resorting  to  higher  precision 
arithmetic.  (We  note  in  passing  that  whereas  Cholesky's  method  produces 
a calculated  LS  solution  at  for  which  the  third  term  in  (D2)  is  of 
negligible  size,  this  is  not  always  the  case  with  some  other  algorithms. 
Hence,  by  ignoring  the  third  term  in  (D2)  these  less  stable  methods  some- 
times introduce  a sizeable  truncation  error  into  (D4).)  Tn  order  to  avoid 
the  use  of  higher  precision  arithmetic,  we  abandon  (D4)  whenever  significance 
losses  become  troublesome  and  resort  instead  to  the  slower  formula  (Dl) . 

When  solving  a sequence  of  problems  for  successive  values  of  m , we  usually 
delay  this  switch  to  subroutine  ACCSSQ  until  the  calculated  values  from 

(D4)  are  about  two  or  three  orders  of  magnitude  less  than  the  quantity 
T 

y y . (Less  experienced  users  may  wish  to  switch  periodically  to  subioutine 
ACCSSQ  in  order  to  gauge  the  accuracy  of  formula  (D4)  for  a particular 


set  of  data.) 


Appo  ini  ix  E ; T(u-  Ko_kn<AN  subprogram 

This  append  1 x contains  listings  ot  six  single  precision  FORTRAN  sub- 
routines which,  togethei  with  an  appropriate  driver  program,  solve  any 
one  of  the  three  lineal  prediction  problems  (1),  (7),  or  (9).  The  user- 

supplied  driver  program  must  define  the  input  parameters  to  subroutine 
. AdNE,  and  arrange  for  display  of  the  output  parameters  from  FABNE.  The 
other  five  subroutines  (CHKAC,  CHSOL,  AKSTOP,  ACCSUM,  and  ACCSS^)  are 
invoked  automat ical ly  by  FABNE,  and  hence  require  no  action  on  the  part 
of  the  user. 

Subroutine  FABNE  forms  the  forward  and/or  backward  prediction  ruirmal 

equations  (6),  (H) , or  (10)  lor  a given  time  series  x, ,x_,...,x  , solve 

— 12  n 

the  normal  equations  by  calls  to  ; ibroutines  t'HFAC  and  CHSOL,  and  deter- 
mines the  "optimal"  number  ot  filter  coefficients  within  the  range 

m sms  m by  calls  to  subroutine  AKSTQP.  Calls  are  also  made 

start  last 

by  FABNE  to  subroutine  ACCSUM  when  forming  inner- products , and  sometimes 
(when  directed  by  the  FABNE  input  parameter  RSCODE)  to  subroutine  ACCSSQ 
when  calculating  the  minimum  residual  sum  of  squares.  FABNE  is  invoked 
by  the  following  calling  sequence : 

CALL  FABNE ( X , N , MSTART , MLAST , ICODE,  AKCODE , RSCODE , MDIM.V, 
VC,S,SC,A,M,SSQ,P,OCODE) 


Description  of  parameters 
Input  parameters 

X - is  a vector  containing  the  time  series. 

N - is  the  length  of  the  time  series. 

MSTART  - is  the  minimum  value  of  m for  which  the  normal 


equations  are  solved. 
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MLAST  - is  the  maximum  value  of  m for  which  the  normal 
equations  could  be  solved. 

ICODE  - is  an  input  code  with  values  as  follows: 

1 - the  forward  prediction  equations  (6)  are  used, 

-1  - the  backward  prediction  equations  (8)  are  used, 

0 - the  forward  and  backward  equations  (10)  are  used. 
AKCODE  - is  an  input  code  with  values  as  follows: 

0 - the  sample  mean  has  not  been  subtracted  out  of  the 

time  series, 

1 - the  sample  mean  has  been  subtracted  out  of  the  time 

series . 

(This  information  is  passed  by  FABNE  to  AKSTOP.) 


RSCODE  - is  an  input  code  with  values  as  follows: 

0 - the  minimum  residual  sum  of  squares  is  calculated 

by  the  fast  formula  (D4) , 

1 - the  minimum  residual  sum  of  squares  is  calculated 

by  subroutine  ACCSSQ,  which  uses  the  slower  (but 
more  accurate)  formula  (Dl) . 

MDIM  - is  used  to  dimension  the  work  space  arrays  and  must  be  set 
to  at  least  MLAST* (MLAST+1 ) /2 . 

Workspace 

V , VC , S , SC  - are  workspace  arrays  used  to  form  (in  V and  S)  and 


solve  (in  VC  and  SC)  the  normal  equations. 


Output  paiametcis 


1H 


A - is  a vectoi  containing  the  "optimal"  filtei  coefficients. 

M - is  the  number  of  "optimal"  filter  coefficients. 

SS\}  - is  the  minimum  residual  sum  of  squares, 
p - is  the  corresponding  mean  square  residual. 

OCODK  - is  an  output  code  with  values  as  follows: 

0 - an  "optimal"  solution  has  been  calculated 

(i.e.  the  usual  exit). 

1 - the  maximum  number  of  filter  coefficients  IMLAST) 

was  used,  but  "optimality"  was  not  achieved. 

2 - premature  teimi nation  has  occurred  because  the 

current  normal  equations  matrix  is  not  |>ositive 
definite  (possibly  due  to  rounding  errors).  In 
this  case  the  previous  solution  is  returned  instead. 

T 

Subroutine  CIIFAO  produces  the  t'holosky  factorization  1.L  of 

any  symmetric  positive  definite  matrix  A . The  lower  triangle  of 

A Is  stored  by  rows  in  the  vector  A , which  is  eventually  overwritten 

T 

by  l,  . Subroutine  CHSOl  solves  Id,  x b , and  hence  it  must  be  pre- 
ceded by  c'HFAC  so  that  the  Cholesky  factorization  Is  already  available. 
The  right-hand  side  b is:  stoied  in  the  vector  R , which  is  eventually 
overwritten  by  the  solution  x . These  two  FORTRAN  subroutines,  which 
are  based  on  the  ALGOL  procedures  choldet.  2 and  cholsol  2 of  Naitin, 
l'oters,  and  Wilkinson  (1971),  were  kindly  supplied  to  us  by  Dr.  M.G.  Cox 
of  the  National  Physical  Laboiatory,  Teddington,  England. 

Subroutine  AKSTOP  implements  the  Akaike  FPE  criterion  described  in 
Appendix  A.  Users  who  wish  to  substitute  some  other  optimality  criterion 


J9 

in  place  of  Akaike's  title  should  have  no  difficulty  in  replacing  AKSTOP 
by  a subroutine  of  their  own  design. 

Subroutine  ACCSUM  calculates  inner-products  accurately,  assuming 
that  rounding  errors  made  in  evaluating  the  products  are  negligible  (see 
Appendix  C) . If  higher  precision  arithmetic  is  allowed,  ACCSUM  could  be 
replaced  by  a higher  precision  subroutine  based  on  the  simpler  program 
segment  (C2)  of  Api>endix  C.  In  most  environments,  savings  in  computer 
time  would  result  from  such  a substitution. 

Subroutine  ACCSS(t  (which  is  only  used  when  RSCODE  ■ 1)  calculates 
the  minimum  residual  sum  of  squares  by  the  slow  formula  (Dl)  of 
Api>endix  D. 

Finally,  a double  precision  version  of  all  six  subroutines  given 
below  can  be  obtained  by  making  the  following  changes: 

(a)  teplace  .ill  REAL  declarations  by  DOUBLE  PRECISION  declarations, 

(b)  teplace  Sv'RT  by  DS^RT  in  lines  17  and  40  of  CHFAC, 

(e)  replace  the  constant  I.E20  by  1.D20  in  the  DATA  statement  in 


line  •>  1 of  FAHNE. 
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SUBROUTINE  FABNE  ( X , N,  MSTART,  MLAST,  ICODE,  AKCODE, 

* RSCODE , Ml)  1 M , V,  VC,  S,  SC,  A,  M,  SSQ,  P,  OCODE ) 

REAL  A,  S,  V,  X,  SC,  VC,  SSQ,  P 

INTEGER  M,  N,  MDIM,  MI.AST,  MSTART,  ICODE,  AKCODE, 

* OCODE , RSCODE 

DIMENSION  X ( N ) , V(MDIM),  VC(MDIM),  S(MLAST), 

* SC  (MI.AST),  A (MI.AST) 

THIS  SUBROUTINE  EORMS  THE  FORWARD  AND/OR  BACKWARD  NORMAL 
EQUATIONS  FOR  AN  N-POINT  TIME  SERIES  X,  SOLVES  THE 
NORMAL  EQUATIONS  BY  CALLS  TO  SUBROUTINES  CHFAC  AND 
CHSOL,  AND  DETERMINES  THE  'OPTIMAL'  NUMBER  OF  FILTER 
COEFFICIENTS  WITHIN  THE  RANGE  MSTART . LE . M . LE . MLAST  BY 
CALLS  TO  SUBROUTINE  AKSTOP.  (IF  THE  REQUIRED  NUMBER  M OF 
FILTER  COEFFICIENTS  IS  KNOWN  BEFOREHAND,  MSTART  AND 
MLAST  SHOULD  BOTH  BE  SET  EQUAL  TO  M . ) THE  NORMAL 
EQUATIONS  ARE  FORMED  EFFICIENTLY  FOR  M-1,2,3,...,  AND 
THEY  ARE  SOLVED  FOR  M “MSTART , MSTART+ 1 , MSTART +2 , . . . UNTIL 
EITHER  M -MLAST  OR  'OPTIMALITY'  IS  ACHIEVED.  CALLS  ARE 
MADE  TO  SUBROUTINE  ACCSUM  WHEN  FORMING  INNER-PRODUCTS, 
AND  SOMETIMES  TO  SUBROUTINE  ACCSSQ  WHEN  CALCULATING  THE 
RESIDUAL  SUM  OF  SQUARES  (SEE  DESCRIPTION  OF  PARAMETER 
RSCODE  BELOW) . 


FORMAL  PARAMETER  LIST: 

INPUT  PARAMETERS 

X A VECTOR  REPRESENTING  THE  TIME  SERIES. 

N THE  DIMENSION  OF  X. 

MSTART  THE  FIRST  VALUE  OF  M FOR  WHICH  THE  NORMAL 
EQUATIONS  ARE  SOLVED  (MSTART . GE . 1 ) . 

MLAST  THE  LARGEST  VALUE  OF  M FOR  WHICH  THE 

NORMAL  EQUATIONS  COULD  BE  SOLVED  (MLAST. GE. 

MSTART);  ALSO  USED  TO  DIMENSION  S , SC , AND  A. 

ICODE  AN  INTEGER  INPUT  CODE  WITH  VALUES: 

-l  - ONLY  BACKWARD  PREDICTION  IS  USED, 
l - ONLY  FORWARD  PREDICTION  IS  USED, 

0 - BOTH  FORWARD  AND  BACKWARD  PREDICTION 
IS  USED. 

AKCODE  AN  INTEGER  INPUT  CODE  USED  IN  SUBROUTINE  AKSTOP 
(SEE  SUBROUTINE  AKSTOP) . 

RSCODE  AN  INTEGER  INPUT  CODE  WITH  VALUES: 

0 - A 'FAST'  CALCULATION  OF  THE  RESIDUAL 

SUM  OF  SQUARES  IS  USED, 

1 - A MORE  ACCURATE  (BUT  SLOWER)  CALCULATION 

OF  THE  RESIDUAL  SUM  OF  SQUARES  IS  USED  BY 
CALLING  SUBROUTINE  ACCSSQ. 

MDIM  THE  DIMENSION  OF  V AND  VC  ( MDIM . GE . MLAST* 
(MLASTFl)/2) . 

WORKSPACE 

V A VECTOR  REPRESENTING  THE  LOWER  TRIANGLE  OF 
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THE  NORMAL  EQUATIONS  MATRIX,  WHICH  IS  STORED 
IN  V BY  ROWS. 

VC  A COPY  OF  V WHICH  IS  USED  IN  THE  CHOLESKY 

FACTORIZATION. 

S A VECTOR  REPRESENTING  THE  RIGHT  HAND  SIDE  OF 

THE  NORMAL  EQUATIONS. 

SC  A COPY  OF  S WHICH  IS  USED  IN  THE  CHOLESKY 

SOLUTION. 

OUTPUT  PARAMETERS 

A A VECTOR  CONTAINING  THE  OPTIMAL  FILTER 

COEFFICIENTS.  \ 

M THE  OPTIMAL  NUMBER  OF  FILTER  COEFFICIENTS. 

SSQ  THE  RESIDUAL  SUM  OF  SQUARES  FOR  THE  OPTIMAL 
M COEFFICIENT  FILTER. 

P THE  MEAN  SQUARE  RESIDUAL  FOR  THE  OPTIMAL 

M COEFFICIENT  FILTER, GIVEN  BY 

P = SSQ  / DENOM, 

2 (N-M)  IF  ICODE.EQ.O 

WHERE  DENOM  * 

(N-M)  IF  ICODE.EQ. (l.OR.-l) . 

OCODE  AN  OUTPUT  CODE  WITH  VALUES: 

0-  ' OPTIMAL ' SOLUTION  FOUND  (I.E.  NORMAL  EXIT), 

1- M. EQ.MLAST,  BUT  'OPTIMALITY'  WAS  NOT 
ACHIEVED  (THE  USER  MAY  WISH  TO  INCREASE 
THE  VALUE  OF  MLAST) . 

2- PREMATURE  TERMINATION,  BECAUSE  THE  CURRENT 
NORMAL  EQUATIONS  MATRIX  IS  NOT  POSITIVE 
DEFINITE  (POSSIBLY  DUE  TO  ROUNDING  ERRORS). 
THE  PREVIOUS  SOLUTION  IS  RETURNED  IN  THIS 
EVENT. 

LOCAL  VARIABLES 

REAL  CHSS , SSQL,  PL,  SUM,  YTY , BIG 

INTEGER  I,  J,  K,  KD , KPM , KPI,  MP1 , NMM , NPl , IFAIL, 

* NMMMl , ITEND 

THE  CONSTANT  BIG  CAN  BE  SET  TO  ANY  LARGE  REAL  NUMBER. 

IT  IS  USED  TO  ASSIGN  AN  INITIAL  VALUE  TO  SSQ  AND  P 
WHEN  MSTART.GT.l. 

DATA  BIG  /1.E20/ 

FORM  AND  SOLVE  THE  NORMAL  EQUATIONS  FOR  M-l . 

OCODE  * 0 
NPl  - N + 1 
M-l 

NMM  -N-M 

CALL  ACCSUM  (X , N,  NMM,  0,  SUM) 


n n o 


42 


IF  (ICODE)  10,  20,  30 
C BACKWARD 

10  YTY  = SUM 

SUM  = SUM  - X ( 1 ) * *2  + X (N ) * *2 
GO  TO  40 

C BOTH 

20  SUM  » 2 . 0 *SUM  - X(l)**2  + X(N)**2 
YTY  * SUM 
GO  TO  40 
C FORWARD 

30  SUM  * SUM 

YTY  = SUM  - X ( 1 ) * * 2 + X(N)**2 
40  V ( 1 ) = SUM 

CALL  ACCSUM (X , N,  NMM , 1,  SUM) 

IF  (ICODE)  50,  60,  70 
C BACKWARD 

50  SUM  = SUM 
GO  TO  80 

C BOTH 

60  SUM  = SUM*2 . 0 
GO  TO  80 
C FORWARD 

70  SUM  =■  SUM 
80  S (1 ) = SUM 

A ( 1 ) = S ( 1 ) /V  ( 1 ) 

SSQL  = YTY  - A ( 1 ) *S  ( 1 ) 

IF  (RSCODE. EQ. 1)  CALL  ACCSSQ(M,  N,  X,  A,  ICODE,  SSQL) 
PL  = SSQL/  (N-M) 

IF  (ICODE. EQ.0)  PL  = SSQL/ ( 2 . 0* (N-M) ) 

IF  (MLAST.EQ.l)  GO  TO  310 
IF  (MSTART.EQ. 1)  GO  TO  90 
SSQL  = BIG 
PL  * BIG 
90  GO  TO  140 

SOLVE  THE  NEW  NORMAL  EQUATIONS  AND  TEST  FOR  OPTIMALITY. 

100  CONTINUE 

KD  = M* (M+l)/2 
IF  (M. LT.MSTART)  GO  TO  140 
KD  = M* (M+l)/2 
DO  110  K*1 , KD 
VC  ( K ) = V (K  ) 

110  CONTINUE 

DO  120  1=1, M 
SC  ( I ) = S ( I ) 

120  CONTINUE 

CALL  CHFAC (M , MDIM , VC,  IFAIL) 

IF  (IFAIL. NE.0)  GO  TO  300 
CALL  CHSOL (M , MDIM,  VC,  SC,  CHSS) 

SSQ  » YTY  - CHSS 

IF  (RSCODE. EQ. 1)  CALL  ACCSSQ(M,  N,  X,  SC,  ICODE,  SSQ) 
P = SSQ/  (N-M) 

IF  (ICODE. EQ.0)  P = SSQ/ (2.0* (N-M) ) 

CALL  AKSTOP  (P , PL,  M,  N,  AKCODE , ITEND) 
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IF  (ITEND.EQ.l)  GO  TO  310 
DO  130  I=1,M 
A(l)  » SC  (I) 

130  CONTINUE 

IF  (M.GE.MLAST)  GO  TO  320 
SSQL  * SSQ 
PL  = P 

FORM  THE  NORMAL  EQUATIONS  FOR  SUCCESSIVE  VALUES  OF  M. 

140  IF  (ICODE)  150,  200,  250 
C BACKWARD 

150  MP1  = M + 1 

K = M * (M*l)/2  ♦ 1 

V (K  ) * S(M)  - X ( 1 ) *X  ( MPl ) 

KPl  = K + 1 

KPM  = K + M 
J = 1 

DO  160  I-KPl,KPM 
J = J + 1 

V(I)  = V(I-MPl)  - X ( J ) *X  (MPl ) 

160  CONTINUE 
K = 0 

DO  180  I = 1 , M 
DO  170  J=1,I 
K = K + 1 

V (K  ) = V (K ) - X (NMM+I ) *X (NMM+J) 

170  CONTINUE 

180  CONTINUE 

DO  190  1=1, M 

S(I)  = S ( I ) - X (NMM)  *X  (NMM  + I ) 

190  CONTINUE 

NMMMl  = NMM  - 1 

CALL  ACCSUM  (X , N.,  NMMMl,  MPl,  SUM) 

S (MPl ) = SUM 

YTY  = YTY  - X (NMM) **2 

M = MPl 

NMM  = N - M 

GO  TO  100 

C BOTH 

200  MPl  ■ M + 1 

K = M* (M+l)/2  + 1 

V (K  ) = S (M)  - X (NMM)  *X  (N  ) - X(1)*X(MP1) 

KPl  = K + 1 

K PM  = K + M 
J = 1 

DO  210  I =KP1 , KPM 
J = J + 1 

V(I)  = V(I-MPl)  - X (NMM)  *X  (NP1-J ) - X(J)*X(MPl) 
210  CONTINUE 
K = 0 

DO  230  1=1, M 
DO  220  J = 1 , 1 

K * K ♦ 1 

V ( K ) * V (K  ) - X(MPl-I)*X  (MP1-J  ) - 

* X (NMM  + I ) *X  (NMM+J) 
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220  CONTINUE 
230  CONTINUE 

DO  240  1=1  , M 

S ( I ) = S ( I ) - Y.  (MPl-I  ) *X  (MP1)  - X (NMM)  *X  (NMM  + I ) 
240  CONTINUE 

NMMM1  = NMM  - 1 

CALL  ACCSUM (X , N,  NMMM1,  MPl , SUM) 

S(MP1)  = 2 . E0  *SUM 

YTY  = YTY  - X (MPl ) * *2  - X(NMM)**2 
M * MPl 
NMM  = N - M 
GO  TO  100 
C FORWARD 

250  MPl  = M + 1 

K = M* (M+l ) /2  + 1 
V (K ) = S (M)  - X (NMM) *X (N) 

KP1  = K + 1 
KPM  = K + M 
J = 1 

DO  260  I =KP 1 , KPM 
J = J + 1 

V ( I ) = V(I-MPl)  - X (NMM)  *X  (NP1-J) 

260  CONTINUE 
K = 0 

DO  280  1 = 1, M 
DO  270  J = 1,I 
K = K + 1 

V(K)  = V (K ) - X (MP1-I) *X (MP1-J) 

270  CONTINUE 
280  CONTINUE 

DO  290  1=1, M 

S ( I ) = S ( I ) - X (MP1-I ) *X (MPl) 

290  CONTINUE 

NMMMl  = NMM  - 1 

CALL  ACCSUM  (X,  N,  NMMMl,  MPl,  SUM) 

S (MPl ) = SUM 
YTY  = YTY  - X (MPl ) * *2 
M = MPl 
NMM  = N - M 
GO  TO  100 

PREPARE  THE  OUTPUT. 

300  OCODE  = 2 

IF  (M.GT.MSTART)  GO  TO  310 
M = 1 

SSQ  = SSQL 
P - PL 
GO  TO  320 
310  SSQ  = SSQL 
P - PL 

IF  (M.EQ.l)  GO  TO  320 
M - M - 1 

320  IF  (ITEND.EQ.O  .AND.  OCODE. NE. 2)  OCODE  - 1 
RETURN 
ND 
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SUBROUTINE  CHFAC (N , NT,  A,  IFAIL) 

REAL  A 

INTEGER  N,  NT,  IFAIL 
DIMENSION  A (NT) 

C 

C CHOLESKY  DECOMPOSITION  LU  (U  - TRANSPOSE  OF  L)  OF 

C N*N  SYMMETRIC  POSITIVE  DEFINITE  MATRIX  A.  LOWER 

C TRIANGLE  OF  A STORED  BY  ROWS  IN  A,  AND  OVERWRITTEN 

C BY  LOWER  TRIANGLE  OF  L.  NT  MUST  BE  NO  SMALLER 

C THAN  N*(N  + l)/2.  BASED  ON  THE  M ART IN -PETERS - 

C WILKINSON  ALGORITHM  CHOLDET2. 

C 

C NO  PROGRAM  UNITS  CALLED. 

C 

C LOCAL  VARIABLES  - 

C 

REAL  ONE,  X,  ZERO,  SORT 
INTEGER  I,  IERROR , J,  K,  P,  0,  R 

CONSTANTS  - 


C 


ZERO  = 0.E0 
ONE  = 1 . OE  0 


IERROR  = 1 
P = 0 

DO  40  I =1  , N 
0 = P + 1 
R = 0 


DO  30  J =1  , I 
X = A ( P + 1 ) 

IF  (Q.GT.P)  GO  TO  20 
DO  10  K =Q , P 


R = R + 1 
X » X - A (K  ) * A (R  ) 

10  CONTINUE 

20  R = R + 1 

P = P + 1 

IF  (I.BQ.J  .AND.  X . LE . ZERO)  GO  TO  50 
IF  (I.EQ.J)  A ( P ) = ONE/SQRT(X) 

IF  (I .NE. J)  A ( P ) = X * A ( R ) 

30  CONTINUE 
40  CONTINUE 

IERROR  * IERROR  - 1 
50  IFAIL  = IERROR 
RETURN 
END 
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SUBROUTINE  CHSOL (N , NT,  L,  B,  SS) 

REAL  B,  L,  SS 
INTEGER  N,  NT 
DIMENSION  B ( N ) , L (NT ) 

SOLVES  LUX  =*  B (U  = TRANSPOSE  OF  THE  N*N  LOWER 
TRIANGULAR  MATRIX  L),  X OVERWRITING  B.  LOWER 
TRIANGLE  OF  L STORED  BY  ROWS  IN  L.  FORMS  IN  SS 
THE  SQUARE  OF  THE  2-NORM  OF  Y - UX.  NT  MUST  BE 
NO  SMALLER  THAN  N*(N  + 1 ) /2 . BASED  ON  THE  MARTIN- 
PETERS-WILKI NS ON  ALGORITHM  CHOLSOL2. 

NO  PROGRAM  UNITS  CALLED 

LOCAL  VARIABLES  - 

REAL  T,  X,  ZERO 

INTEGER  I,  IREV,  K,  KREV,  P,  0.  S 
CONSTANTS  - 
ZERO  = 0.E0 

SOLVE  LY  * B AND  FORM  SS  - 

T = ZERO 
P = 1 

DO  30  I = 1 , N 
0 = 1-1 
X = B ( I ) 

I F (0- LT. 1)  GO  TO  20 
DO  10  K=1 ,0 

X = X - L(P)*B(K) 

P « P + 1 
10  CONTINUE 

20  B ( I ) = X*L(P) 

P » P + 1 
T = T + B ( I ) * *2 
30  CONTINUE 
SS  = T 

SOLVE  UX  = Y - 

DO  60  IREV-1  ,N 

I « N + 1 - IREV 
S ■ P - 1 
P « S 
0 = 1 + 1 
X = B ( I ) 

IF  (Q.GT.N)  GO  TO  50 
DO  40  KREV»Q,N 

K - N + 0 - KREV 
X - X - L (S ) * B (K  ) 

S - S - K + 1 
40  CONTINUE 

50  B ( I ) = X *L  (S ) 

60  CONTINUE 
RETURN 
END 


a. 
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SUBROUTINE  AKSTOP ( P , PL,  M,  N,  AKCODE , ITEND) 

REAL  P,  PL 

INTEGER  M,  N,  AKCODE,  ITEND 
C 

C FOLLOWING  THE  AKA  1 KE  GOODNESS-OF -F IT  CRITERION,  WE 
C SAY  THAT  M-l  IS  OPTIMAL  IF  M IS  THE  SMALLEST  POSITIVE 
C INTEGER  FOR  WHICH, 

C 

C IF  THE  MEAN  IS  SUBTRACTED  OUT, 

C 

C P* (N+ (M+l ) )/(N- (M+l ) ) .GT.PL*  (N+-M)  / (N-M) , 

C 

C BUT  OTHERWISE, 

C 

C P* (N+M) /(N-M) .GT.PL* (N+ (M-l) ) /(N- (M-l) ) . 

C 

C 

C HERE  N IS  THE  LENGTH  OF  THE  TIME  SERIES, 

C M IS  THE  NUMBER  OF  COEFFICIENTS,  P IS  THE  MEAN  SQUARE 
C RESIDUAL,  AND  PL  IS  THE  CORRESPONDING  QUANTITY  FOR 
C M-l  COEFFICIENTS. 

C 

C AKCODE  IS  AN  INTEGER  INPUT  PARAMETER  WHICH 
C MUST  BE  SET  TO  1 I F THE  MEAN  IS  SUBTRACTED 
C OUT  OF  THE  TIME  SERIES,  AND  WHICH  MUST  BE 
C SET  TO  0 OTHERWISE. 

C ITEND  IS  AN  INTEGER  OUTPUT  PARAMETER  WHICH  IS  SET 
C TO  1 IF  THE  AKA  IKE  CRITERION  IS  SATISIFIED,  AND  WHICH 
C IS  SET  TO  0 OTHERWISE. 

C 

C LOCAL  VARIABLES  - 
REAL  XN,  XM 
ITEND  = 0 
XN  = N 
XM  * M 

IF  ( AKCODE. NE. 1)  GO  TO  10 

IF  (P* ( (XN+XM+1 . 0) /(XN-XM-1 . 0) ) .GT.PL* ( (XN+XM) / (XN-XM) 

* ) ) ITEND  = 1 
RETURN 

10  CONTINUE 

IF  (P*( (XN+XM) / (XN-XM) ) .GT.PL* ( ( XN+XM-1 . 0 ) / ( XN-XM+1 . 0 ) 

* ) ) ITEND  = 1 
RETURN 

END 
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SUBROUTINE  ACCSUM(X,  N,  NMM,  MOVE,  SUM) 

REAL  X,  SUM 
INTEGER  N,  NMM,  MOVE 
DIMENSION  X (N  ) 

C 

C CALCULATES  INNER-PRODUCTS  ACCURATELY,  ASSUMING  THAT 
C ROUNDING  ERRORS  MADE  IN  EVALUATING  THE  PRODUCTS  ARE 
C NEGLIGABLE. 

C 

C LOCAL  VARIABLES  - 
C 

REAL  S,  C,  Y,  T 
INTEGER  I 
S =>  0 .0 
C = 0 .0 
DO  10  1*1, NMM 

Y = C + X (I)*X (I+MOVE) 

T = S + Y 
C = (S-T)  + Y 
S = T 
10  CONTINUE 
SUM  = S + C 
RETURN 
END 
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SUBROUTINE  ACCSSQ(M,  N,  X,  A,  ICODE,  SSQ) 
REAL  X,  A,  SSQ 
INTEGER  M,  N,  ICODE 
DIMENSION  X (N)  , A(M) 

CALCULATES  SUM  OF  SQUARES  ACCURATELY,  FROM  FIRST 
PRINCIPLES. 

LOCAL  VARIABLES  - 
REAL  APX , TSSQ 
INTEGER  I,  J,  NMM 
TSSQ  = 0.0 
NMM  = N - M 
IF  (ICODE)  50,  10,  10 
C FORWARD  ONLY 
10  CONTINUE 
SSQ  = 0.0 
DO  30  1=1, NMM 
APX  = 0.0 
DO  20  J = 1 , M 

APX  = APX  + A(J) *X  (M+I-J) 

20  CONTINUE 

SSQ  = SSQ  + (X (M+I )-APX) **2 
30  CONTINUE 

IF  (ICODE)  50,  40,  80 

C SAVE  FORWARD  SSQ  FOR  'BOTH'  CALCULATION. 

40  TSSQ  = SSQ 
C BACKWARD  ONLY 
50  CONTINUE 
SSQ  =0.0 
DO  70  1=1, NMM 
APX  =0.0 
DO  60  J=1  ,M 

APX  = APX  + A(J  ) *X  (J+I  ) 

60  CONTINUE 

SSQ  = SSQ  + (X (I )-APX) **2 
70  CONTINUE 

C ADD  FORWARD  AND  BACKWARD  TO  FORM  'BOTH'  SSQ. 

SSQ  = SSQ  + TSSQ 
80  CONTINUE 
RETURN 
END 


